# Note: To execute this code outside the 02_Replicate.R, load packages by
# uncommenting #library() commands, set working directory to replication folder 
# using setwd(). Then load the datasets by uncommenting the commands loading the 
# data.

#library("dplyr")
#library("ggplot2")
#library("lfe")
#library("estimatr")
#library("stargazer")
#library("mvtnorm")

#setwd() # set working directory to replication folder
#load("Data/rtl_data.Rdata")
#lh <- read.csv("Data/laver_hunt_party_org.csv")
#params <- read.csv(file = "Data/params.csv")
#output <- read.csv(file = "Data/simulation_output.csv")


# Table 1 -----------------------------------------------------------------

m1 <- felm(lead_center ~ losepower + ternary_class_factor  + voteshare + outofcoalition|
             0|0|party, data = dat, weights = dat$weights1[dat$cs == 1], subset = cs == 1)
m2 <- felm(lead_center ~ losepower + ternary_class_factor  + voteshare + outofcoalition|
             elect|0|party, data = dat, weights = dat$weights1[dat$cs == 1], subset = cs == 1)
m3 <- felm(lead_center ~ losepower + ternary_class_factor  + voteshare + outofcoalition|
             party|0|party, data = dat, weights = dat$weights1[dat$cs == 1], subset = cs == 1)
m4 <- felm(lead_center ~ losepower + ternary_class_factor  + voteshare + outofcoalition|
             decade + party|0|party, data = dat, weights = dat$weights1[dat$cs == 1], subset = cs == 1)
m5 <- felm(lead_center ~ losepower + ternary_class_factor  + voteshare + outofcoalition|
             party+elect|0|party, data = dat, weights = dat$weights1[dat$cs == 1], subset = cs == 1)

# Regression table
stargazer(m1, m2, m3, m4, m5, 
          keep = "losepower",
          dep.var.labels = "Center Platform$\\_{t+1}$", 
          covariate.labels = "Loss of Power$\\_t$", 
          dep.var.caption = "",
          add.lines = list(c("Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes"),
                           # c("$\\Delta$ Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes"),
                           c("Platform$\\_t$ FE", "yes", "yes", "yes", "yes", "yes"),
                           c("Party FE", "", "", "yes", "yes", "yes"),
                           c("Decade FE", "", "", "", "yes", ""),
                           c("Election FE", "", "yes", "", "", "yes")),
          omit.stat = c("f", 
                        "ser",
                        "rsq",
                        "adj.rsq"),
          float = F,
          out = "Outputs/Table_1.tex")

# Table 2 -----------------------------------------------------------------

n0 <-  felm(to_extreme ~ losepower + to_extreme_t_1  + voteshare + outofcoalition|
              0|0|party, data = dat,weights = dat$weights4)
n1 <- felm(to_extreme ~ losepower * to_extreme_t_1  + voteshare + outofcoalition|
             0|0|party, data = dat,weights = dat$weights4)
n2 <- felm(to_extreme ~ losepower * to_extreme_t_1+ ternary_class_factor + voteshare + outofcoalition|
             elect|0|party, data = dat,weights = dat$weights4)
n3 <- felm(to_extreme ~  losepower * to_extreme_t_1 + ternary_class_factor  + voteshare + outofcoalition|
             party|0|party, data = dat,weights = dat$weights4)
n4 <- felm(to_extreme ~  losepower * to_extreme_t_1 + ternary_class_factor  + voteshare + outofcoalition|
             decade + party|0|party, data = dat,weights = dat$weights4)
n5 <- felm(to_extreme ~  losepower * to_extreme_t_1 + ternary_class_factor  + voteshare + outofcoalition|
             party+elect|0|party, data = dat,weights = dat$weights4)



stargazer(n0, n1, n2, n3, n4, n5, 
          keep = c("losepower", "to_extreme_t_1", "losepower:to_extreme_t_1"),
          dep.var.labels = "Center Platform$\\_{t+1}$", 
          covariate.labels = c("Loss of Power$\\_t$", "To Extreme$\\_{t-1}$", "Loss of power$\\_t \\times$ To Extreme$\\_{t-1}$"), 
          dep.var.caption = "",
          add.lines = list(c("Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes", "yes"),
                           #      c("$\\Delta$ Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes", "yes"),
                           c("Platform$_t$ FE", "yes", "yes", "yes", "yes", "yes", "yes"),
                           c("Party FE", "", "", "", "yes", "yes", "yes"),
                           c("Decade FE", "", "", "", "", "yes", ""),
                           c("Election FE", "", "", "yes", "", "", "yes")),
          omit.stat = c("f", 
                        "ser",
                        "rsq",
                        "adj.rsq"),
          float = F,
          out = "Outputs/Table_2.tex")

# Table 3 -----------------------------------------------------------------

# The two panels of Table 3 are combined manually. This code creates the tables from
# which Table 3 was constructed.

o1 <- felm(absolute_shift ~ losepower * ls + absolute_shift_t_1 + 
             voteshare + outofcoalition|0|0|party, data = dat, weights = dat$weights5)

o2 <- felm(absolute_shift ~ losepower * ls + absolute_shift_t_1 + 
             voteshare + outofcoalition|elect|0|party, data = dat, weights = dat$weights5)

o3 <- felm(absolute_shift ~ losepower * ls + absolute_shift_t_1 + 
             voteshare + outofcoalition|party|0|party, data = dat, weights = dat$weights5)

o4 <- felm(absolute_shift ~ losepower * ls + absolute_shift_t_1 + 
             voteshare + outofcoalition|party + decade|0|party, data = dat, weights = dat$weights5)

o5 <- felm(absolute_shift ~ losepower * ls + absolute_shift_t_1 + 
             voteshare + outofcoalition|party + elect|0|party, data = dat, weights = dat$weights5)


q1 <- felm(to_extreme ~ losepower * ls + to_extreme_t_1 + 
             voteshare + outofcoalition|0|0|party, data = dat, weights = dat$weights5)

q2 <- felm(to_extreme  ~ losepower * ls + to_extreme_t_1 + 
             voteshare + outofcoalition|elect|0|party, data = dat, weights = dat$weights5)

q3 <- felm(to_extreme  ~ losepower * ls + to_extreme_t_1 + 
             voteshare + outofcoalition|party|0|party, data = dat, weights = dat$weights5)

q4 <- felm(to_extreme  ~ losepower * ls + to_extreme_t_1 + 
             voteshare + outofcoalition|party + decade|0|party, data = dat, weights = dat$weights5)

q5 <- felm(to_extreme  ~ losepower * ls + to_extreme_t_1 + 
             voteshare + outofcoalition|party + elect|0|party, data = dat, weights = dat$weights5)



stargazer(o1, o2, o3, o4, o5, 
          keep = c("losepower", "ls", "losepower:ls"),
          dep.var.labels = "Shift Magnitude$\\_{t+1}$", 
          covariate.labels = c("Loss of Power$\\_t$", "Large Selectorate", "Loss of power$\\_t \\times$ Large Selectorate$"), 
          dep.var.caption = "",
          add.lines = list(c("Voteshare$\\_t$", "yes",  "yes", "yes", "yes", "yes"),
                           #     c("$\\Delta$ Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes"),
                           c("Out of Coalition$\\_{t}$", "yes", "yes", "yes", "yes", "yes"),
                           c("S$_t$ FE", "yes", "yes", "yes", "yes", "yes"),
                           c("Party FE",  "", "", "yes", "yes", "yes"),
                           c("Decade FE",  "", "", "", "yes", ""),
                           c("Election FE",  "", "yes", "", "", "yes")),
          omit.stat = c("f", 
                        "ser",
                        "rsq",
                        "adj.rsq"),
          float = F,
          out = "Outputs/Table_3A.tex")



stargazer(q1, q2, q3, q4, q5, 
          keep = c("losepower", "ls", "losepower:ls"),
          dep.var.labels = "To Extreme$\\_{t+1}$", 
          covariate.labels = c("Loss of Power$\\_t$", "Large Selectorate", "Loss of power$\\_t \\times$ Large Selectorate$"), 
          dep.var.caption = "",
          add.lines = list(c("Voteshare$\\_t$", "yes",  "yes", "yes", "yes", "yes"),
                           c("Out of Coalition$\\_{t}$", "yes", "yes", "yes", "yes", "yes"),
                           c("S$_t$ FE", "yes", "yes", "yes", "yes", "yes"),
                           c("Party FE",  "", "", "yes", "yes", "yes"),
                           c("Decade FE",  "", "", "", "yes", ""),
                           c("Election FE",  "", "yes", "", "", "yes")), 
          omit.stat = c("f", 
                        "ser",
                        "rsq",
                        "adj.rsq"),
          float = F,
          out = "Outputs/Table_3B.tex")


# Table 4 -----------------------------------------------------------------


m1 <-  felm(lead_govt_party2 ~ losepower * TE + voteshare + outofcoalition|
              elect|0|party, subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])
m2 <-  felm(lead_govt_party2 ~ losepower * TE +voteshare + outofcoalition|
              party|0|party, subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])
m3 <-  felm(lead_govt_party2 ~ losepower * TE + voteshare + outofcoalition|
              elect + party|0|party, subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])
m4 <-  felm(lead_govt_party2 ~ losepower * AS + voteshare + outofcoalition|
              elect|0|party,  subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])
m5 <-  felm(lead_govt_party2 ~ losepower * AS + voteshare + outofcoalition|
              party|0|party, subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])
m6 <-  felm(lead_govt_party2 ~ losepower * AS + voteshare + outofcoalition|
              elect + party|0|party, subset = cs == 1, data = dat, 
            weights = dat$weights6[dat$cs == 1])



stargazer(m1, m2, m3, m4, m5, m6,
          keep = c("losepower", "TE", "losepower:TE", "AS", "losepower:AS"),
          dep.var.labels = "Government Party$\\_{t+2}$",
          covariate.labels = c("Loss of Power$\\_t$", 
                               "To Extreme$\\_{t+1}$", 
                               "Absolute Shift$\\{t+1}$",
                               "Loss of power$\\_t \\times$ To Extreme$\\_{t+1}$",
                               "Loss of power$\\_t \\times$ Absolute Shift$\\_{t+1}$"), 
          dep.var.caption = "",
          add.lines = list(c("Voteshare$\\_t$", "yes", "yes", "yes", "yes", "yes", "yes"),
                           c("Platform$_t$ FE", "yes", "yes", "yes", "yes", "yes", "yes"),
                           c("Party FE", "", "yes", "yes", "", "yes", "yes"),
                           c("Election FE", "yes", "", "yes", "yes", "", "yes")),
          omit.stat = c("f", 
                        "ser",
                        "rsq",
                        "adj.rsq"),
          float = F,
          out = "Outputs/Table_4.tex")

# Figure 1 ----------------------------------------------------------------

# Subset to parties prior to election t
ingovt <- filter(dat, main_sample == 1)


l_post <- lm_robust(I(lead_ternary_class == -1) ~ losepower, data = ingovt, clusters = party)
l_pre <- lm_robust(I(ternary_class == -1) ~ losepower, data = ingovt, clusters = party)
m_post <- lm_robust(I(lead_ternary_class == -0) ~ losepower, data = ingovt, clusters = party)
m_pre <- lm_robust(I(ternary_class == -0) ~ losepower, data = ingovt, clusters = party)
r_post <- lm_robust(I(lead_ternary_class == 1) ~ losepower, data = ingovt, clusters = party)
r_pre <- lm_robust(I(ternary_class == 1) ~ losepower, data = ingovt, clusters = party)

# This function extracts relevant estimates from the regression models above
get_cond_means <- function(model, platform, t){
  means = coef(model)[1] + c(0, coef(model)[2])
  ses = c(summary(model)$coef[1,2],
          sqrt(vcov(model)[1,1] + vcov(model)[2,2] + 2 * vcov(model)[1,2]))
  z = c("Remained in power", "Lost Power")
  return(data.frame(platform, t, means, ses, z))
}

ests <- do.call("rbind",
                list(
                  get_cond_means(l_post, "Left platform", 1),
                  get_cond_means(l_pre, "Left platform", 0),
                  get_cond_means(m_post, "Center platform", 1),
                  get_cond_means(m_pre, "Center platform", 0),
                  get_cond_means(r_post, "Right platform", 1),
                  get_cond_means(r_pre, "Right platform", 0))) %>%
  as.data.frame() %>%
  mutate(Platform = factor(platform, levels = unique(platform))) 
br =  c(1:5/10)
labs = c(paste0(1:5*10, "%"))

fig1 <- 
  ggplot(ests, aes(x = t, y= means, col = z)) +
    geom_point(position = position_dodge(width =0.15)) +
    facet_grid(~Platform) +
    geom_line( position = position_dodge(width =0.15)) +
    geom_errorbar(aes(ymin = means - 1.96 *ses, ymax = means + 1.96 * ses), width = .1,
                         position = position_dodge(width =0.15)) +
    geom_errorbar(aes(ymin = means - 1.64 *ses, ymax = means + 1.64 * ses), 
                         position = position_dodge(width =0.15), lwd = 1.1, width = 0) +
    scale_x_continuous("Election", breaks = c(0,1), labels = c("Current", "Next")) + 
    theme_minimal() + 
    scale_color_manual("Result of current election", values = c("#2c7fb8", "#7fcdbb")) + 
    scale_y_continuous("Proportion of parties adopting...", breaks =br, labels = labs) +
    theme(legend.position = "bottom")

ggsave(fig1, file = "Outputs/Figure_1.pdf", width = 9, height = 5)



# Figure 2 ----------------------------------------------------------------

ingovt <- filter(dat, main_sample == 1)

m1 <- lm_robust(lead_govt_party2 ~ losepower * TE, data = dat, clusters = party)
m2 <- lm_robust(lead_govt_party2 ~ losepower * TE, data = ingovt, clusters= party)

get_cis <- function(model){
  betas <- rmvnorm(n = 10000, mean = coef(model), sigma = vcov(model))
  x_mat1 <- cbind(1, seq(-0.83, 0.88, by = .001))
  x_mat2 <- cbind(1, 1,  seq(-0.83, 0.88, by = .001), seq(-0.83, 0.88, by = .001))
  nolp <- betas[,c(1,3)]
  
  
  nolp_pred <-x_mat1 %*% coef(model)[c(1,3)]
  lp_pred <- x_mat2%*% coef(model)
  nolp_ci <- t(apply(nolp %*% t(x_mat1), 2, FUN = quantile, probs = c(0.025, 0.975)))
  lp_ci <- t(apply(betas %*% t(x_mat2), 2, FUN = quantile, probs = c(0.025, 0.975)))
  return(data.frame(TE = rep(seq(-0.83, 0.88, by = .001),2),
                    losepower = rep(c("Did not lose power", "Lost power"), each = length(nolp_pred)),
                    pred = c(nolp_pred,lp_pred), 
                    ci_lo = c(nolp_ci[,1], lp_ci[,1]),
                    ci_hi = c(nolp_ci[,2], lp_ci[,2])))
}

t <- get_cis(m1) %>% mutate(sample = "All parties")
t2 <- get_cis(m2) %>% mutate(sample = "Governing parties prior to election t")

rug_dat <- data.frame(TE = c(dat$TE, ingovt$TE),
                      sample = c(rep("All parties", nrow(dat)),
                                 rep("Governing parties prior to election t", nrow(ingovt))))
labs <- data.frame(TE = c(-.38, .4),
                   pred = c(.95, .95),
                   labs = c(sprintf('\u2190 To center'), "To extreme \u2192"))


fig2 <- 
  rbind(t, t2) %>%
    filter(TE  > -0.5 & TE < 0.5) %>% 
    ggplot(aes(x = TE)) +
    facet_grid(~sample) +
    geom_line(lwd = 1.1, aes(y = pred, col = losepower)) +
    theme_minimal() +
    geom_ribbon(aes(ymin = ci_lo, ymax = ci_hi, fill = losepower), alpha = 0.4) +
    scale_y_continuous("Probability of being governing party\nafter election t+1") +
    scale_x_continuous("Shift to extreme from election t to election t+1", limits = c(-.5,.5),
                       breaks = c(-2:2/4), labels = c(-50, -25, 0, 25, 50)) +
    geom_rug(data = rug_dat, sides = "b") + 
    scale_color_manual("Result of election t", values = rev(c("#2c7fb8", "#7fcdbb"))) +
    scale_fill_manual("Result of election t", values = rev(c("#2c7fb8", "#7fcdbb"))) +
    geom_vline(xintercept = 0, lty = 2) + 
    geom_text(data = labs, aes(y = pred, label = labs), cex = 3) +
    theme(legend.position = "bottom")

ggsave(fig2, file = "Outputs/Figure_2.pdf", width = 9, height = 5)


# Figure 3 ----------------------------------------------------------------

# Note that Figure 3 is an illustration and does not use data.


# Figure 4-6 ----------------------------------------------------------------

# These figures are produced from the Markov chain simulation output. To run
# the simulation, please use the Mathematica file 05_simulation.nb.

# construct a dataframe with exogenous parameters and Markov chain output
names(output) <- c("RRwext", "RRlext", "RLwext", "RLlext", "LRwext", "LRlext", "LLlext", "LLwext",
                   "RPrRev", "LPrRev", "RRincWadj", "RRincLadj", "RLincWadj", "RLincLadj",
                   "LRincWadj", "LRincLadj", "LLincWadj", "LLincLadj", "Rprobvic")
names(params) <- c("xm", "sigma", "gammav")

mat <- cbind(params, output)


fig4 <- ggplot(mat, aes(x = xm, y = sigma, fill = RRlext - RRwext, col = RRlext - RRwext)) +
  facet_grid(.~gammav, labeller = label_bquote(cols = lambda == .(gammav))) + geom_raster() +
  scale_fill_gradientn(expression(paste("Difference in \nPr(Extreme) following \nLoss vs. Victory")),
                       colours = c("red","white","blue"),
                       breaks = c(-.02, 0, .02), labels = c(-.02, 0, .02)) + 
  geom_vline(xintercept = 0) + theme_minimal() +
  scale_x_continuous(expression(paste(y[m], " (Location of Median Voter)"))) +
  scale_y_continuous(expression(paste(pi, " (Selectorate Size)")))
ggsave(fig4, file = "Outputs/Figure_4.pdf", width = 8, height = 3)



fig5 <- ggplot(mat, aes(x = xm, y = sigma, fill =RPrRev)) +
  facet_grid(.~gammav, labeller = label_bquote(cols = lambda == .(gammav))) + geom_raster() +
  scale_fill_gradientn(expression(paste("Difference in \nPr(Reversal) following \nLoss vs. Victory")),
                       colours = c("red","white","blue"),
                       breaks = c(-.015, 0, .015), labels = c(-.015, 0, .015)) + 
  geom_vline(xintercept = 0) + theme_minimal() +
  scale_x_continuous(expression(paste(y[m], " (Location of Median Voter)"))) +
  scale_y_continuous(expression(paste(pi, " (Selectorate Size)")))

ggsave(fig5, file = "Outputs/Figure_5.pdf", width = 8, height = 3)


fig6 <- ggplot(mat, aes(x = xm, y = sigma, fill =RRincLadj - RRincWadj)) +
  facet_grid(.~gammav, labeller = label_bquote(cols = lambda == .(gammav))) + geom_raster() +
  scale_fill_gradientn(expression(paste("Difference in \nPr(Adjustment) following \nLoss vs. Victory")),
                       colours = c("red","white","blue"),
                       breaks = c(-.05, 0, .05), labels = c(-.05, 0, .05)) + 
  geom_vline(xintercept = 0) + theme_minimal() +
  scale_x_continuous(expression(paste(y[m], " (Location of Median Voter)"))) +
  scale_y_continuous(expression(paste(pi, " (Selectorate Size)")))

ggsave(fig6, file = "Outputs/Figure_6.pdf", width = 8, height = 3)
